Solid molecular hydrogen: The Broken Symmetry Phase 
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By performing constant-pressure variable-cell ab initio molecular dynamics simulations we find 
a quadrupolar orthorhombic structure, of Pca2\ symmetry, for the broken symmetry phase (phase 
II) of solid H2 at T = and P = 110 — 150 GPa. We present results for the equation of state, 
lattice parameters and vibronic frequencies, in very good agreement with experimental observations. 
Anharmonic quantum corrections to the vibrational frequencies are estimated using available data 
on H2 and D2. We assign the observed modes to specific symmetry representations. 
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The quest for the structure of the high-pressure phases 
of hydrogen is a long standing one. Early predictions of 
an insulator- metal transition ffl, led to a large body of 
experimental and theoretical work during the past sixty 
years. Metallization was not found as promptly as ini- 
tially expected pj but a rich phase diagram emerged. The 
current picture of the low and room-temperature phase 
diagram up to ~ 230 GPa is essentially based on the 
optical studies performed in diamond anvil cell (DAC) 
devices during the past decade HQ . 

There is a consensus that hydrogen exhibits at least 
three different solid molecular phases. 1) At low pres- 
sures (< 110 GPa for para-H^) the centers of the H2 
molecules crystallize into an hep lattice, but zero-point 
motion overcomes rotational energy barriers, leading to 
a free-rotator phase (phase I). 2) Between 110 and 150 
GPa intermolecular interactions freeze the molecular ro- 
tations into an ordered broken-symmetry phase (BSP, or 
phase II). 3) Above ~ 150 GPa, a third phase (H-A, or 
phase III) is attained, whose structure is unknown. Here 
we focus specifically on phase II. Optical measurements 
in phase II indicate the presence of two, possibly three, 
infrared (IR) active modes, in constrast to phase I, where 
only one mode is observed. Raman spectra show a sin- 
gle peak, at a frequency ~ 10 % lower than that of the 
IR modes. A consistent picture of the structural and 
dynamical properties of phase II is still lacking || . 

On the theoretical side, most of the existing work con- 
sists of static total energy calculations within the local 
density approximation (LDA). Zero-point energy (ZPE) 
of the protons has been in a few cases included a posteri- 
ori based on frozen phonon calculations . Hexagonal 
close packed structures with two and four molecules per 
unit cell appear to be the strongest candidates for the 
ground structure of phase II, but the relative orientation 
of the molecules is uncertain [^|-[Io|. Cubic structures 
were also suggested [[|. This uncertainty persists, due to 
the incomplete optimization of the lattice parameters and 
the a priori selection of a particular space group symme- 
try in static calculations. Ab initio Molecular Dynamics 
simulations, which do not rely on the choice of a specific 
space group, have also been reported |ll|,Q. However, 
these earlier attempts were hampered by a fixed simula- 
tion cell and a poor Brillouin Zone (BZ) sampling. 

We report extensive ab initio Molecular Dynamics sim- 
ulations of solid H2 in the 110 — 150 GPa range. Addi- 
tional crucial ingredients are (i) a variable cell, constant 
pressure approach that largely eliminates prejudice on 
the space group symmetry, molecular orientations and 
lattice constants jl3]; (ii) fully converged BZ sampling, 
achieved by a freshly implemented k ■ p technique p[ . 
Electronic correlations are treated within the LDA sup- 
plemented with Becke-Perdew gradient corrections [ p"5| , 
while the proton-electron interaction is described through 
a local pseudopotential Jl6[ | requiring an energy cutoff of 
60 Ry fl7|] . We stress that electrostatic and band energies 



are fully included in this approach |l(| . 

Preliminary runs using large simulation cells (up to 
128 atoms) with T-point sampling only, produced layered 
ground state structures with in-plane triangular ordering 
of molecules. However, we found a strong dependence of 
the results on cell shape and atom number, confirming 
that an accurate BZ sampling is crucial p0 18 1. 
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vergence on BZ sums for insulating molecular H2 turned 
out to require at least 8x8x8 fc-points in the full BZ of 
a four-molecule cell 1 10 . This is particulary demanding 



in our approach because of the lack of symmetry con- 
straints. Thus, at variance with the standard approach 
where the electronic orbitals for each fc-point in the BZ 
are expanded on a generic basis set (e.g. plane waves), 
we adopt a k ■ p expansion in terms of the (occupied plus 
empty) self-consistently generated T-point orbitals: 
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where uf(r) is the periodic part of the i-th Bloch function 
at wave number k, and a^. are unitary matrices obtained 
by direct diagonalization the k ■ p hamiltonian at k [ fl9| . 
This expansion is strictly equivalent to the standard ap- 
proach when the number of states ./V equals the size of 
the basis set [|o|. In practice the sum can be truncated 
at a significantly smaller number of T-point orbitals with 
excellent accuracy [Q. With this improved functional 
we studied a system of 8 atoms with 512 fc-points in the 
full BZ plj ], and with N = 64. We also performed sim- 
ulations with 32 atoms, 128 fc-points and N = 96, and 
with 64 atoms, 64 fc-points and N = 128. 

The 8-atom simulation started from the Pmc2\ struc- 
ture proposed by Kaxiras and Broughton Q (see Fig. 
1), that is the most favorable among a class of structures 
containing also the P2/m. Pressure was set to 140 GPa. 
When constant-pressure dynamics was switched on, the 
system spontaneously transformed into Kitaigorodskii 
and Mirskaya's quadrupolar structure Pca2i |22| , pro- 
posed for H 2 by Nagara and Nakamura || (see Fig. 1). 
We then performed additional runs starting with differ- 
ent initial structures, namely the cubic Pa3 and the or- 
thorhombic Cmca (see Fig. 1), previously proposed for 
solid H2 at lower and higher pressures, respectively |Pl[l8|- 
Although these structures can be be ruled out for phase 
II on the basis of the number of vibronic modes j^] , they 
represent reasonable starting points for our dynamical 
structural search. In the range of existence of phase II, 
the Cmca structure turned out to be unstable against 
Pca2i, while Pa3 was found to be locally stable, but 
much higher in enthalpy than Pcet2 1 . The quadrupo- 
lar structure of symmetry P2\jc j^] was also consid- 
ered and found to be locally stable, with an enthalpy 
0.015 eV/molecule higher than that of the Pca2i (at 140 
GPa). The simulations were repeated with supercells of 
32 atoms (double size along the two in-plane directions), 
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and 64 atoms (double size in all directions). However, we 
failed to observe any signature of ground state structures 
with unit cells larger than 8 atoms, and Pca2y prevailed 
in all cases. 

The above results have been obtained by treating the 
protons as classical particles, since quantum effects have 
been proven irrelevant to the relative hierarchy of molec- 
ular phases j23j. However, given the high level of accu- 
racy of the present calculations, it seems wise to include 
the quantum effects. By rotating a single molecule on the 
Pca2\ structure p4| we do not find multiple energy wells, 
suggesting that tunnelling (of high order in K) should 
not be important, at least near 150 GPa. Quantum cor- 
rections were therefore evaluated a posteriori within the 
harmonic approximation (i.e. to first order in Ti). To 
this aim, we generated molecular dynamics trajectories 
at very low temperature (T « 2 K), from which we ex- 
tracted [^5| the frequencies of the 21 modes (4 vibrons, 
9 phonons and 8 librons) corresponding to the T-point 
of the 4-molecule BZ. The ZPE, computed by summing 
over the above modes, was evaluated only for structures 
Pca2i and P2y/c, the closest competitors for phase II. 
In fact, a Pa3 space group is incompatible with exper- 
imental data 0, while all remaining structures are un- 
stable in the clamped-nuclei approximation. The insta- 
bility implies the existence of a distortion with negative 
clamped-nuclei energy curvature, and excludes stabiliza- 
tion through quantum corrections, to first order in h. 
The ZPE turned out to be lower in P2i/c by only 0.003 
eV/molecule with respect to Pca2y, yielding a quantum- 
corrected enthalpy difference of 0.012 eV/molecule, still 
in favor of the classicaly determined Pca2y structure. 

It is worth stressing that in Pca2y the centers of the H2 
molecules occupy the sites of a slightly distorted hep lat- 
tice. The temperature-induced transformation of phase 
II into phase I should then be accompanied by orienta- 
tional disordering and by full recovery of the hep symme- 
try. Thus, the equation of state (EOS) of phases I and II 
are expected not to be significantly different. In Fig. 2 we 
report our EOS and the c/a ratio calculated for the Pca2i 
structure in the range from 100 to 150 GPa. The b/a ra- 
tio turns out to be essentially pressure-independent, al- 
though slightly smaller (b/a « 1.715) than the hep value 
of y/3. Comparison with the experimental results for the 
EOS and the c/a ratio in phase I [^6| (full line) is also 
rather satisfactory. 

A structure with Pca1\ symmetry is compatible with 
optical measurements on phase II, where a strong Ra- 
man peak, and two IR modes are observed §,§. A sym- 
metry analysis of the four vibronic modes of the Pca2y 
structure, usually labeled as Ay, A 2 , By, and B 2 [p|,p7|, 
shows that all of them should be Raman active, and 
three of them (Ay, By, and B 2 ) IR active. The frequen- 
cies of these four modes were evaluated by decomposing 
molecular dynamics trajectories at 100 and 150 GPa into 
the symmetry representations of space group Pca2y p8[. 



The resulting frequencies turned out to be substantially 
higher than experiment. Moreover, they increased as a 
function of pressure, in agreement with the calculated re- 
duction of the molecular bond length (see Table I), but 
in flagrant disagreement with the observed trends, par- 
ticularly for the Raman mode. 

A closer glance at the experimental data, however, in- 
dicates that quantum effects are far from negligible (see 
Fig. 7 of Ref. ||). In fact, the D2 Raman frequency 
(scaled by \/2) is roughly 200 cmT 1 higher than that of 
H2. To first order in h, quantum corrections to the fre- 
quency are known to scale as the inverse of the particle's 
mass, no matter what the potential is pjjfl . Thus, using 
the experimental H 2 and D 2 frequencies, we obtain em- 
pirical but accurate quantum corrections to the classical 
values following Ref. [3C[| . In the case of the H 2 Raman 
mode the correction is as large as 500 cm -1 , while for 
the IR modes the effect is more modest, amounting to 
about 150 cm -1 . The resulting quantum-corrected vi- 
bron frequencies are reported in Fig. 3, and compared 
with experimental data on H 2 ||. The agreement with 
experiment is now very good, and shows that the Raman 
vibron turnover is a purely quantum effect. Moreover, 
the agreement allows us to interpret the experimental 
Raman mode as our Ay mode: a totally symmetric, in- 
phase vibration of the four molecules in the unit cell. Ac- 
cording to its symmetry classification, this mode should 
also be IR active. It does not appear in the experimental 
IR spectra, possibly due to a small oscillator strength. 
This is reasonable for a mode where all molecules oscil- 
late in phase, since IR activity can only come from a 
small intramolecular asymmetry. 

We attribute the two IR vibrons to the By and B 2 
modes. The predicted splitting of the two IR vibrons 
turns out to be in fairly good agreement with the data. 
Finally, we predict a Raman vibron, of symmetry A 2 , 
with a "classical" frequency similar to that of the IR 
modes, being similarly an out-of-phase oscillation of the 
four molecules. Applying the same type of quantum cor- 
rection as for the IR modes, the A 2 should in principle be 
observed in the region around 4500 cm -1 . Failure to ob- 
serve modes A 2 , By and B 2 in Raman spectroscopy could 
be attributed to their small Raman strength, contribu- 
tions from out of phase vibrations canceling out. The 
observed Raman mode could alternatively be attributed 
to mode A 2 , but this would imply a quantum-corrected 
frequency of 4370 cm" 1 (at 150 GPa), in clear disagree- 
ment with experiment. Symmetry lowering arguments 
and larger unit cells have been very recently invoked to 
explain a larger number of optical modes by band folding 
in ortho-para D 2 mixtures |l|. The ortho -para distinc- 
tion is however not included in our calculation. 

Finally, our calculation provides the full electronic 
structure and pressure coefficients, and in principle the 
optical properties. In Fig. 4 we report the band structure 
and the electronic density of states of the Pca2i struc- 
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ture, at P = 100 GPa (the density of states at 150 GPa 
is also reported, for comparison). A gap of about 4 eV 
opens up the otherwise fairly free-electron- like bands. As 
shown (see also Table I), the gap decreases monotonically 
with pressure but is always finite and large in the range 
of stability of phase II. 

In summary, by means of extensive and accurate 
constant-pressure ab-initio molecular-dynamics simula- 
tions we have obtained what appears to be a definitive 
description of the broken symmetry phase of H2. The 
Pca2\ structure, first described by Kitaigorodskii and 
Mirskaya ]23] for the ground state of a classical quadrupo- 
lar system, and later proposed for H2 ||, is favored. 
Equation of state and vibronic frequencies turn out to 
be in very good agreement with available experimental 
data. 

This work was partly sponsored by INFM project LO- 
TUS, and by the European Commission under contract 
ERBCHRXCT 940438. 
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FIG. 1. Candidates for the ground state structure of solid 
H2 in the broken symmetry phase (phase II), projected along 
the c-axis. Black (gray) arrows represent molecules centered 
on the c (c/2) plane and pointing towards the positive-.? hemi- 
sphere. 

FIG. 2. Calculated equation of state of H2 in phase II: 
8-atom simulation cell (solid circles), and 64-atom simulation 
cell (open squares). The experimental EOS of H2 in phase 
I [26] is also reported (solid line) . The inset shows the pressure 
dependence of the c/a ratio (triangles are data from [26], other 
symbols as above). 

FIG. 3. IR and Raman vibron frequencies of H2 as a func- 
tion of pressure. The lower and the couple of upper solid lines 
are the experimental Raman and IR data [3,5], respectively. 
Circles are calculated Raman vibrons (mode Ai). Squares and 
hexagons are calculated IR active vibrons (modes B\ and B2 , 
respectively). Open triangles are estimates for the new A2, 
Raman-only vibron. 

FIG. 4. Electronic band structure (left panel) and density 
of states (right panel) of solid H2 in the Pca2\ structure, at 
P = 100 GPa. The energy gap is 4.12 eV. The density of 
states at 150 GPa is also reported (dashed line, right panel) 
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TABLE I. Bond length, volume, and band gap in the two 
quadrupolar structures Pca2\ and P2i/c, as a function of 
pressure. 





P (GPa) 


d H H (A) 


V (cm 3 /mole) 


E g (eV) 


Pca2 1 


100 


1.377 


2.75 


4.12 


Pca2 1 


150 


1.373 


2.35 


2.53 


P2 1 /c 


150 


1.374 


2.35 


2.61 



Pca2 1 P^/c 



Pmc^ Cmca 

t 4 t 
4 t 4 
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Pressure (GPa) 



